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ABSTRACT 

Aims. We investigate the behaviour of nonlinear, nonideal Alfven wave propagation within an inhomogeneous magnetic environment. 
Methods. The governing MHD equations are solved in ID and 2D using both analytical techniques and numerical simulations. 
Results. We find clear evidence for the ponderomotive effect and visco-resistive heating. The ponderomotive effect generates a 
longitudinal component to the transverse Alfven wave, with a frequency twice that of the driving frequency. Analytical work shows 
the addition of resistive heating. This leads to a substantial increase in the local temperature and thus gas pressure of the plasma, 
resulting in material being pushed along the magnetic field. In 2D, our system exhibits phase mixing and we observe an evolution in 
the location of the maximum heating, i.e. we find a drifting of the heating layer 

Conclusions. Considering Alfven wave propagation in 2D with an inhomogeneous density gradient, we find that the equilibrium 
density profile is significantly modified by both the flow of density due to visco-resistive heating and the nonlinear response to the 
localised heating through phase mixing. 
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' 1. Introduction 

. Phase mixing, a mechanism for dissipating shear Alfven waves, 
' was first proposed by Heyvaerts & Priest (1983). The basic con- 
cept is straightforward: consider shear Alfven waves propagat- 
ing in an inhomogeneous plasma, such that on each magnetic 
fieldline each wave propagates with its own local Alfven speed. 
] Thus, after propagating a certain distance, these neighbouring 
perturbations will be out of phase, which will lead to the gener- 
ation of smaller and smaller transverse spatial scales, and thus 
, to the growth of strong cuiTents. This ultimately results in a 
■ strong enhancement of the dissipation of Alfven-wave energy 
via viscosity and/or resistivity. Alfven wave phase mixing has 
also been studied extensively as a possible mechanism for heat- 
ing the corona (e.g. Heyvaerts & Priest 1983, Browning 1991, 
Ireland 1996, Malara et al. 1996, Narain & Ulmschneider 1990; 
1996). 

Nakariakov et al. (1997) extended the model of Heyvaerts & 
Priest (1983) to include compressibility and nonlinear effects. 
They found that fast magnetoacoustic waves are generated con- 
tinuously by Alfven-wave phase mixing at a frequency twice that 
of the driven Alfven wave, and that these generated waves prop- 
agate across magnetic fieldlines and away from the phase mix- 
ing layer. Since there is a permanent leakage of energy away 
from the phase mixing layer, these fast waves can cause indi- 
rect heating of the plasma as they propagate away and dissipate 
far from the layer itself, thus spreading the heating across the 
domain. Nakariakov et al. (1997) found that the amplitude of 
these fast waves grows linearly in time, according to weakly- 
nonlinear, j6 = 0, analytical theory. Nakariakov et al. (1998) fur- 
ther extended this model to include a background steady flow 
and found that the findings of their 1997 model persist. 
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Hood et al. (2002) investigated the phase mixing of single 
Alfven-wave pulses and found that this results in a slower power- 
law damping (as opposed to the standard ~ exp(-f^) for har- 
monic Alfven waves, i.e. Heyvaerts & Priest 1983). Hence, 
Alfvenic-pulse perturbations will be able to transport energy to a 
greater coronal height than that of a harmonic Alfven wavetrain. 
Botha et al. (2000) considered a developed stage of Alfven-wave 
phase mixing and found that the growth of the generated fast 
waves saturates at amplitudes much lower than that of the driven 
Alfven wave. They concluded that the nonlinear generation of 
fast waves (Nakariakov et al. 1997) saturates due to destruc- 
tive wave interference, and has little effect on the standard phase 
mixing model of Heyvaerts & Priest (1983). The numerical sim- 
ulations of Botha et al. assumed wave propagation in an ideal 
plasma, and thus heating was absent from their model. Tsiklauri 
and co-authors repeated these numerical experiments for both a 
weakly and strongly-nonlinear Alfvenic pulse (Tsiklauri et al. 
2001; 2002). 

De Moortel et al. (1999; 2000) investigated how gravitational 
density stratification and magnetic field divergence changes the 
efficiency of phase mixing. They report that the resultant dissipa- 
tion can be either enhanced or diminished depending on the spe- 
cific choice of equilibrium. Ruderman et al. (1998) investigated 
phase mixing in open magnetic equilibria, and found similar re- 
sults using the WKB method. These 1998 analytical results were 
repeated and corrected by numerical and semi-analytical work 
by Smith et al. (2007). All these authors found the following: a 
diverging magnetic field enhances the efficiency of phase mix- 
ing, whereas gravitational stratification diminishes the mecha- 
nism. Hood et al. (1997a; 1997b) derive analytical, self-similar 
solutions of Alfven-wave phase mixing in both open and closed 
magnetic topologies. 

Ofman & Davila (1995) found that in an inhomogeneous coro- 
nal hole with an enhanced dissipation parameter (S = 1000 - 
10,000), Alfven waves can dissipate within several solar radii, 
which can provide significant energy for the heating and acceler- 
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ation of the solar wind. This model was later extended to include 
nonlinear effects (Ofman & Davila 1997). 
Parker (1991) pointed out that phase mixing requires an ignor- 
able coordinate, an assumption which is expected to be unphys- 
ical in the corona. Parker found that including all three co- 
ordinates results in the driven Alfven waves coupling with a 
fast magnetoacoustic mode, and that this elimates the growth 
in transverse spatial scales, and thus phase mixing is absent 
from the system. Instead, such a system exhibits resonant ab- 
sorption (e.g. Lee & Roberts 1986; HoUweg & Yang 1988) on 
the surfaces where the phase velocity equals the Alfven veloc- 
ity. However, Parker's conclusions are in disagreement with the 
work of Tsiklauri and co-authors who considered the interac- 
tion of an impulsively-generated, weakly-nonlinear MHD pulse 
with a one-dimensional density inhomogeniety, considered in 
the three-dimensional regime (i.e. without an ignorable coordi- 
nate) in both an ideal (Tsiklauri & Nakariakov 2002) and re- 
sistive (Tsiklauri et al. 2003) plasma. Tsiklauri and co-authors 
found that phase mixing remains a relevant paradigm and that 
the dynamics can still be quahtatively understood in terms of 
the classic 2.5D models. Mocanu et al. (2008) have revisited 
the Heyvaerts & Priest (1983) model using anisotopic viscosity 
(i.e. incorporating the Braginskii 1965 stress tensor) and report 
that this significantly increases the damping lengths, i.e. com- 
pared to those obtained for isotropic dissipation. More recently, 
Threlfall et al. (2010) have investigated the effect of the Hall 
term on phase mixing in the ion-cyclotron range of frequencies. 
Another key concept that we shall invoke in this paper is the 
ponderomotive force: a nonhnear force proportional to spatial 
gradients in magnetic pressure, also referred to as the Alfven 
wave-pressure force. The ponderomotive effect has been con- 
sidered in a solar context initially by Hollweg (1971) and later 
by Verwichte and co-authors (Verwichte 1999; Verwichte et 
al. 1999). Hollweg (1971) considered linearly-polarised Alfven 
waves propagating in a direction parallel to the magnetic field, 
and found that the transverse behaviour of the Alfven wave was 
identical when comparing the linear and nonlinear (to second 
order) solutions, but that longitudinal wave velocity and density 
fluctuations appear in the nonlinear solutions driven by gradi- 
ents in the wave magnetic-field pressure, i.e. the Alfven wave 
is no longer purely transverse and is compressive (through non- 
hnear couphng to magnetoacoustic waves). Thus, the pondero- 
motive force can be used as an extra acceleration mechanism 
and as an explanation for density fluctuations in the solar wind. 
Verwichte (1999) presented a mechanical analogy for the pon- 
deromotive effect by considering the resulting motion of discrete 
particles (beads) on an oscillating string. Verwichte et al. (1999) 
considered the temporal evolution of a weakly-nonhnear, Alfven 
wave in a = homogeneous plasma. These authors showed 
that the an initially-excited, gaussian-pulse perturbation in trans- 
verse velocity splits into two Alfven wave pulses, each propagat- 
ing in opposite directions (as naturally expected). Furthermore, 
Verwichte et al. found that the ponderomotive force produces a 
shock in longitudinal velocity at the starting position. Note that 
in a cold plasma, there is no force to counteract this ponderomo- 
tive acceleration. 

In this paper, we investigate the nonlinear, nonideal behaviour 

of Alfven-wave propagation and phase mixing over long 
timescales. Thus, this paper can be seen as an extension of the 
model of Botha et al. (2000) to include visco-resistive effects. 
We also seek to address a fundamental question of phase mix- 
ing: by considering nonlinear, nonideal phase mixing over long 
timescales, is it possible to observe a drifting of the heating 
layer? In other words, phase mixing will occur due to the den- 



sity inhomogenity in our system, and this process will gener- 
ate strong, localised heating due to enhanced dissipation. This 
locahsed heat deposition is expected to modify the equilibrium 
density profile, and thus may change the location of maximum 
heating. However, it is unclear what the actual result will be: we 
may observe a change in the location of the heating layer, or the 
phase mixing mechanism may break down, or the heating may 
bifurcate spatially (since our density profile will no longer be 
monotonic). Is is also unclear how this may affect the indirect 
heating of the plasma, due to the coupling to the fast magne- 
toacoustic mode (Nakariakov et al. 1997). Such an investigation 
requires nonhnear and nonideal effects to be considered and, as 
we shall see, observed over long timescales. 
The work of Ofman et al. (1998) is also relevant here. These 
authors investigated a model of resonant absorption that incor- 
porated the dependence of loop density on the heating rate, and 
studied the spatial and temporal dependence of the heating layer. 
Ofman et al. find that the heating occurs in multiple resonance 
layers, rather than the single layer of the classic resonant ab- 
sorption models (e.g. lonson 1978; Ulmschneider et al. 1991; 
Ruderman & Roberts 2002) and that these layers drift through- 
out the loop to heat the entire volume. Poedts & Boynton (1996) 
also investigated resonant absorption using nonlinear, resistive 
MHD simulations and found a spreading of the heat deposition, 
i.e. a broadening of the resonant layer due to changes in the back- 
ground inhomogeneity. 

This paper has the following outline: §2 describes the governing 
equations, assumptions and analytical and numerical details of 
our investigation, §3 investigates the ID nonhnear, nonideal sys- 
tem (with no density inhomogeneity) and focuses on the under- 
lying physical processes. §4 details the bulk-flow phenomenon 
found to be present in our system, and investigates its density 
dependence. §5 considers a 2D model with an inhomogeneous 
density profile, and details the long-term evolution and coupled 
nature of the three MHD waves present in our system. The con- 
clusions are presented in §6 and there are two appendicies. 



2. Basic Equations 

We consider the nonlinear, compressible, viscous and resistive 
MHD equations appropriate to the solar corona: 



^+V.(pv) = 0, 
ot 



dt 



H- (V ■ V) V 



' V X B\ 
-Vp + I — I X B -I- vV • S , 



— = V X (v X B) -H TyV^B 
ot 



(1) 



dp 

-^ + iv-V)p = -ypV ■y+ 

ot a 



^ ^Ljl' + v(r-i)e. 



where p is the mass density, v is the plasma velocity, B the 
magnetic induction (usually called the magnetic field), p is the 
plasma pressure, ju = 4;r X 10"^Hm"' is the magnetic permeabil- 
ity, cr is the electrical conductivity, 77 = l/yucr is the magnetic 
diffusivity, y = 5/3 is the ratio of specific heats and j = V x B/// 
is the electric current density. In this paper, vj and v are assumed 
to be constants. 

The viscous stress tensor, S, is given by 
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where e is the rate-of-strain tensor, given by: 

1 / dvi dvi 



2 \dxj dxij 
The viscous heating term is given by 



Qv 



3 3 



1=1 j=i 



Note that in the presence of strong magnetic fields, the classi- 
cal viscous term used in equations (1) is not the most appro- 
priate for the solar corona, since the viscosity takes the form 
of a non-isotropic tensor. The mathematical details of this non- 
isotropic tensor can be found in Braginskii (1965). In this paper, 
we have chosen to implement the simplest, isotropic version of 
the stress tensor as a first step in investigating nonlinear phase 
mixing, and non-isotropic viscosity will be considered in future 
investigations. Our choice of viscosity term also allows a direct 
comparison with the results of Heyvaerts & Priest (1983) and al- 
lows us to derive analytical solutions to our governing equations. 
We now consider a change of scale to non-dimensionalise all 
variables. Let v = vqv*, B = BB*, x = Lx*, y = Ly*, p = pop*, 
P = Pop*, V = = TAt*, 7] - rjo and v = vq, where * denotes 

a dimensionless quantity and vq, B, L, po, po, ta, iJo and vq are 
constants with the dimensions of the variable they are scahng. 
We then set B/ yfjlpo - vq and vq - L/ta (this sets vq as a 
constant equilibrium Alfven speed). We also set tjqTa/L? = R^, 
where is the magnetic Reynolds number. This process non- 
dimensionalises equations (1) and under these scalings, f* = 1 
(for example) refers to f = ta = L/vo', i-e. the time taken to 
travel a distance L at the background Alfven speed. For the rest 
of this paper, we drop the star indices; the fact that all variables 
are now non-dimensionahsed is understood. 



2.1. Basic equilibrium 

We consider equations (1) in Cartesian coordinates and assume 
there are no variations in the z-direction (d/dz = 0). We con- 
sider an inhomogeneous plasma in the x-direction embedded in 
a uniform magnetic field in the y-direction, i.e.: 

Po = P()(x) , Po = constant , Bo = (0, Bo, 0) , vq = (2) 

with finite amplitude perturbations of the form: 



p = p(x,y,t) , v = (v,r,Vy,v,) = v(x,>',f) 
p = p(x,y,t) , B = (B^,By,B,) = B(x,y,t) 

2.2. Analytical description 



(3) 



Following the analysis of Nakariakov et al. (1995; 1997) and 
Botha et al. (2000), equations (1) are written in component form 
using equilibrium (2) and perturbations (3): 





dB^ dy^ 
-df-^'^ 


-Ri 


= Ns 




dBy ^ ^ dY^ 
dt dx 


-Ri 


= Ne 




dB, dw, 
-df-^'-d^ 




= Nj 


dp 
'dt 




dVy\ 

dyj 


= Ns 



(8) 

(9) 
(10) 

(11) 



where the left-hand-sides involve only linear terms and Ni-N% 
denote the nonlinear terms. The nonlinear terms are: 



N2 



N3 = - 



A'4 



Bz dB, 

IX dx ^ dt 

By I dBy dB; 

fx \ dx dy 

B, dB, dvy 
fi dy ^ dt 
B^ I dBy dB,, 
IX \ dx dy 
dv, , J d 
'dt 

By dB, B^ dB^ 
II dy n dx ' 

d 



dx 



(po+p)|v.-+v,-^v. 



-p 



(P0+P)^V.-+V,-^V. 



— {y^By - YyB^) , 



N5 

Ne = --^[vxBy-VyB^) , 

A^8 



dx 
d_ 
dx 



d d 

(V,B, - V,,B,) + — {w,By - VyB,) 



dv^ <9v. 



(12) 

(13) 

(14) 

(15) 
(16) 
(17) 
(18) 
(19) 



d_ d_ 
'dx ^dyj"^ "^\dx dy 

the resistive terms (where Ri-R^ are hnear, R^ is nonUnear) are: 



Ri 
Ri 
R3 



d^ 52 \ 



idx^^df] 



B, 



R4 ^ iy-D- 



dB, 



dBy dBA^ IdB,'^ 
dx dy j \dx 







dp d ^ , d^y 
^ + ^(pov.)+po^ 


= Ni , 


(4) 


Vi = 


(9v.i 
'dt 


dp 
dx 


^BoldBy_dBA^^^ 
H \dx dy 1 


= N2 , 


(5) 


V2 = 






dVy dp 


= N3 , 


(6) 


V3 = 






dvz Bo dB^ 

Po-^ 5~ + ^3 

dt p. dy 


= N4, 


(7) 


V4 = 



and the viscous terms (Vi-Vs are linear, V4 is nonlinear) are 

(2 d^v, dh^ 1 d\, 
-v\---^ + 



3 dx^ 

y 



dy^ 

dh, 2 d^Vy 



3 dxdy ) 
1 d^vA 



dx^ 3 dy"^ 3 dxdy j 
cJ^v, cJ^v 



' dx^ dy'^ 



= v(r-i) 




(20) 
(21) 
(22) 
(23) 

(24) 
(25) 
(26) 



1 Idvx _^ dvy^ 
2\dy dx 
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1 /5v, 
"^2 i dx 



dyj 2\dy 



dx 



dy 



dx dy 



are quadratic in terms of B, and so will have an amplitude that is 
related to the square of the amplitude of the linear Alfven wave. 

(27) 2.3. Numerical Set-up 



Equations (4 - 11) describe all the waves occurring in the plasma. 
These equations can be combined to obtain the evolution expres- 
sions for the velocity components: 







1 


dVi 


PO 








^^ 


^ dx 



dx^ 

_m 

dx 

dR2 
dx 



,2 " 

^dy^ 
dNi 

_ dN5 



2 d\' 
dxdy 



dRi 

'dy 



d^ -,d^\ 2 '^^v.t 



,dfi 



•dy 

Po \ dt 

ydfi ^^dy^ 
1 



dVj 

'~dr 



dxdy 
dNj dNs 
~dt dy 



dN4 

IT 



IdNi 


dR3\ 




"^1 



(28) 



(29) 



(30) 



where v\ix) = B^/fipo(x) is the unperturbed Alfven speed and 
Cj(x) = JPo/poix) is the unperturbed sound speed. In equations 
(28 - 30) the nonlinear and nonideal terms can be found on the 
right-hand-side. Note that, equations (4 - 30) reduce to the cor- 
responding equations in Botha et al. (2000) in an ideal plasma. 
Considering equation (30) and neglecting the nonlinear and non- 
ideal terms (right-hand-side), we see that v, corresponds to the 
linear Alfven wave. Thus, if v^ = Va(x), then both v, and B. will 
depend on x. Hence, there will be coupling to both magnetoa- 
coustic modes as follows: Firstly, we can see that the first term 
in N2 (i.e. ~ B^dB-Jdx) will generate a fast wave (i.e. drives v^^ 
term) via phase mixing. Secondly, we see that the first term in A^3 
(i.e. ~ B.dB,/dy) corresponds to the Alfven- wave pressure gra- 
dient, or the ponderomotive force, and that this (nonlinear) force 
generates a slow wave (i.e. drives v^ term) along the background 
magnetic field. Also note that these two terms we have described 



We solve MHD equations (1) using the LARE2D Lagrangian- 
remap code (Arber et al. 2001). 

We consider a numerical domain with -100 < x < 100, 
< y < 10,000 with uniform-grid spacing in the jc-direction 
and a stretched grid in the y-direction. The stretched grid places 
the majority of the y-direction gridpoints near low values of y. 
Results presented in this paper have a typical numerical resolu- 
tion of 2,000 X 20,000 which means that 6x x 6y ~ O.l (due 
to the stretched grid). We chose such a large domain to ensure 
that no wave energy is reflected back into our system during our 
numerical simulations. 

We drive our system with linearly-polarised, sinusoidal Alfven 
waves. This harmonic wave train is driven at the y = boundary 
such that: 



v^{x, 0, f) = A sin (ojt) , v^ix, 0, f) - \y(x, Q,t) -Q 



(31) 



B^ix, 0, f) = -A ^jpix) sin (cat) , — — 



dB, 

'dy 



dBy 

'dy 







where A is the amplitude and to is the frequency. All other 
quantities have zero gradient boundary conditions at the 
y-boundaries, and we utilise periodic boundary conditions at the 
jc-boundaries. In this paper, we set A - 0.01 and oj = 0.25. At 
the start of the numerical run, the driven Alfven wave is ramped 
up to A over the first four wavelengths. It should be noted that 
boundary conditions (31) do not drive a pure Alfven wave into 
the numerical domain; slow magnetoacoustic modes are nonlin- 
early excited as well with amplitude of order A- (see §3.1). 
Note that our choice of equilibrium (equation 2) is formally only 
an equilibrium in ideal MHD. However, numerical tests with 
no driving term show that this equilibrium is still valid over 
long timescales, i.e. the results that follow are due to the driven 
Alfven waves, not due to our choice of equilibrium. 

3. One-dimensional system 

Let us start by considering the nonlinear and nonideal aspects 
of the Alfven wave. In this section, we set po - constant, i.e. 
we remove the density inhomogeneity. This effectively reduces 
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Fig. 2. (a) Transverse perturbations V; for rj = 0.01, y6 = 0.1 plasma (v = 0), where dashed lines represent the envelope V; = ±Aexp(A;/y). (b) 
Longitudinal perturbations v,, for same plasma, where dashed line represents y = c j and blue horizontal line denotes v,, = 0. In both subfigures, 
t = IOOOt^ and dotted line represents y = Va?. 



equations (4 - 30) to a one-dimensional (ID) system where vari- 
ations in the jc-direction are ignorable (d/dx = 0). This ID sys- 
tem will clearly illustrate the nonlinear driver in the equations as 
well as the effect of viscosity and resistivity. These effects must 
be identified before trying to interpret the results with a density 
inhomogeneity. 

The driven, linear Alfven wave has the form: 

V, = T(aj[t - y/^A]) , B,^- ^/JI^v, , (32) 

where the arbitrary function ;F can include a form for the ramp- 
up of the driver as well as the periodic driver. To illustrate the 
terms, we ignore the ramp-up and consider the sinusodial driver; 
the solutions are: 

Vj. - A sin ((L)[t - y/\ a]) for < y < Vyif , 

= - y/Jip^A sin (w[f - y/v^]) for < y < v^f , 

and 

Vj - B, = for y > v^f . 

The nonlinear evolution of the Alfven wave in an ideal, /3 - 
plasma can be found in Appendix A. However, this result can be 
easily derived from the /3 case (§3.1) and so Appendix A is 
only included for completeness. 

3.1. Ideal /3 > plasma 

Consider a /? plasma, where we set j8 = 0.1. The addition of 
finite-/? has no effect on the v- component of the driven Alfven 
wave (i.e. equation 32 is still valid) but does influence the longi- 
tudinal component of velocity generated by the ponderomotive 
force. This can be seen in Figure la. Here, we see that there are 
two kinds of longitudinal motions in the system: the first can be 
seen between y = cj ^ 280 and y - v^f - 1000. These lon- 
gitudinal motions have a maximum amplitude of (v^ - Cj) 
(this maximum is denoted by the horizontal blue line) and they 
are always positive. Thus, there is a net outflow. 
The second type of longitudinal motion can be seen between 
y = and y - c j. These longitudinal perturbations are acous- 
tic waves (in ID) with both positive and negative motions of 
the plasma along the field, and can be identified as a boundary- 
driven slow wave. Note that the longitudinal motions for < y < 



Cgt are a combination of these two motions: the acoustic / slow 
wave (speed Cs) and the ponderomotive perturbations (speed v^). 
In Figure lb, we can see a blow-up of the region 100 < x < 200 
for Vy (black line). The red line represents an analytical solution 
describing both types of longitudinal motions (see equations 33 
below). The agreement is excellent. Note that the slow wave has 
only propagated a distance of y = Cjf at f = IOOOta and so, after 
y - Cst ~ 280, only one wave is present. 

The analytical solution for these nonlinearly-generated longi- 
tudinal motions is found by substituting equation (32) into the 
weakly-nonlinear form of equation (29), namely: 

dfi dy^ ) ~ po dt ~ jupo dtdy \ 2 / ' 
Thus, the analytical solution for v,, is: 

-^4^ (cos[2w(f - ylc,)] 
_J -cos[2w(f-y/vA)]) 

5^1^ (I - cos[2w(f - y/vA)]) c,t<y< v^f , 

y > Va? . 

which is valid only for early times (i.e. until nonlinear response 
becomes non-negligible). Note that these longitudinal motions 
were first reported in Botha et al. (2000). 

Note that for < y < c^t, both types of wave in equation (33) 
must have the same frequency (i.e. twice the driving frequency) 
but will propagate at different speeds. Thus, they must have dif- 
ferent wavenumbers and hence different wavelengths. For the 
parameters chosen in this paper, these wavelengths are: 

^alfven -12.6 and /Islow = 3.63 . 

Thus, for the typical numerical resolutions considered in this pa- 
per (Sy ~ 0.1), both these perturbations are well resolved. 

3.2. Resistive plasma 

Let us now consider a resistive plasma, where we set rj = 0.01 
(and V -Q,/5 - O.I). As before, an Alfven wave is driven into the 
numerical domain and the resultant v^ behaviour can be seen in 
Figure 2a. It clear that the wave now experiences resistive damp- 
ing. Such damping can be estimated from Fourier analysing the 
linear form of equation (30) to give the dispersion relation: 

= (v^ -H//7w)fc2 (34) 
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Fig. 3. (a) Longitudinal perturbations v, for r; = 0.01, y6 = 0.1 plasma comparing = 0.01, v = (blue) and v = 0.01, rj = Q (red). Note only 
< y < 350 is shown and that the agreement after y » 280 is excellent (the two curves lie on each other), {b) Blow-up of v,, over region Q <y < 100 
for V = 0.01, T] = Q plasma only, where crosses indicate grid point resolution. Black line indicates Vj. = 0. In both subfigures, t = lOOOr^i. 




where w = 0.25, y\ - |Bo| /po, ?? - 0.01 and k = kR + ik/. 
In Figure 2a, the dashed lines represent the envelope Vj. - 
+A exp (kjy), where kj < 0. Assuming » rjai (or equivalently 
that the magnetic Reynolds number is large) we can estimate 
that: 



(35) 



to a high degree of accuracy. 

Figure 2b shows the corresponding longitudinal motions (v,) in 
our resistive system. As in Figure lb, we see that the acoustic 
and ponderomotive wave components are again present, but that 
there is now a third phenomenon: a bulk flow in the positive 
y-direction with v,, > and that has a maximum around y = c j. 
Physically, this bulk flow (i.e. movement of material) is due to 
the ohmic heating in our system, and is a natural consequence of 
driving an Alfven wave in a resistive plasma. 
The longitudinal behaviour is governed by equation (29), which 
in ID (i.e. d/dx - 0) and at early times (i.e. no nonlinear feed- 
back) reduces to: 



yupo at \ dy 



(y-l)?7 d 
/jpo dy 



dy 



(36) 



The right-hand-side of this equation has two contributions: the 
first term depends upon the rate-of-change of the magnetic pres- 
sure gradient ■^,{~f) of the driven wave and the second depends 
upon pressure gradients created by the ohmic heating term (i.e. 

increasing the gas pressure. Note the second contribution 
vanishes under the ideal approximation. 

Let us now evaluate these terms under the assumption that our 
linear Alfven wave can be represented as: 



Vj. = A sin {(jjt - ki(y)e'"-' . 

Thus, from equations (7) and (10), B~ has the form: 

(jL)B(,A r, . , , 

v2(^2 + k]) 

-ki cos {cDt - kRy)] e'"^ . 



(37) 



(38) 



The derivation of equation (38) is given in Appendix B. 
Hence, the right-hand-side of equation (36) can be simplified 
such that: 
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where the terms representing the Ponderomotive force, P, the 
resistive heating force, R, and the non-oscillatory resistive bulk 
flow force, BF, are given by: 



+kism2(ojt-kRy)] e^'^ 
R = -^!^^^^[kRsin2icot-kRy) 

+ki cos2 (cot -kuy)] e^'"^ 
o/A^ (y - 1) r]ki ^2k,y 



BF 



(39) 

(40) 
(41) 



which is valid for c^t < y < V/if. The derivation of these equa- 
tions is given in Appendix B. 

Note that the amplitude of each term is proportional to the am- 
plitude squared of the driven Alfven wave and that they decay 
with height since ki < 0. In addition, we note that P and R are 
trigonometric terms with twice the driven frequency and twice 
the wavenumber. In the absence of dissipation (77 = and so 
ki = 0) there is no decay and R and BF vanish. 
We also note that BF is, critically, non-oscillatory and it is this 
term that will create the bulk flow in the longitudinal direction. 

3.3. Viscous plasma 

Let us now consider a purely viscous plasma, where we set v = 
0.01, T] - and yS — 0.1. Again, we drive an Alfven wave into 
our numerical domain and the resultant v. is identical to that 
seen Figure 2a. However, the damping mechanism is now due 
to viscosity rather than resistivity, as in §3.2. Such damping can 
be estimated from Fourier analysing the linear form of equation 
(30) to give the dispersion relation: 



- + ivoj^ k^ 



(42) 



Note the similarity to equation (34). 

Figure 3a shows a comparison of Vy for the purely-resistive 
(blues, §3.2, Figure 2b) and purely-viscous (red line) systems. 
Focusing on the viscous propagation, we see that three compo- 
nents are present: the ponderomotive wave (between 280 <y < 
1000), the bulk flow in the positive 3'-direction, and the acous- 
tic perturbations (from Q < y < 280). However, note that the 
acoustic perturbations now take a different form to those in the 
purely-resistive system (black line) and we note that in the vis- 
cous system the smaller wavelengths are rapidly damped out by 
viscosity. Above y » 280, the agreement is exceUent (the two 
curves lie ontop of one another). Note that the oscillation is fuUy 
resolved, as can be seen in Figure 3b, which shows a blow-up of 
Vy in the region < >■ < 100 in the viscous system. 
In the viscous plasma, the bulk-flow phenomenon now comes 
from viscous heating, as opposed to ohmic heating as in §3.2. 
Comparing the two systems (Figure 3 a) we see that the agree- 
ment for the ponderomotive component and bulk flow is excel- 
lent (unsurprising, since equations 34 and 42 are identical for 
V = rj). However, the viscous damping of the acoustic compo- 
nent is more pronounced (with an estimated damping length of 
13.5). 

4. Bulk-Flow Phenomenon 

Let us now investigate the effect that our longitudinal bulk flow 
has on the equilibrium density profile, starting with the ideal 



case. Figure 4a shows the evolution of the density perturbations 
&Xt = IOOOta in our ideal plasma. We can clearly see that the 
contributions from the ponderomotive component and acoustic 
wave components. We also see that the density profile is per- 
turbed about the equilibrium but that these perturbations do not 
grow in time. 

Figure 4b shows the density profile at f = lOOOr^i in our visco- 
resistive {tj - v - 0.01) system. Here, in addition to the contri- 
butions from the ponderomotive component and acoustic wave, 
we see there is also a third component; i.e. the bulk flow has 
shifted the density profile in the positive y-direction, leading to 
a decrease in the density profile around y - 0. 

The physical interpretation of this bulk flow is as follows: visco- 
resistive heating increases the local temperature (with maximum 
around y - 0). This increase in temperature increases the ther- 
mal pressure, and the resultant pressure gradient drives the bulk 
flow. Note that the bulk flow is not due to the ponderomotive 
wave pressure force, otherwise it would be apparent in our ideal 
plasma (i.e. in Figure 4a). Thus, the bulk flow is a natural conse- 
quence of driving an Alfven wave in a nonideal plasma. 
Let us now look look at the long-term evolution of the density 
profile. Figure 5a shows the evolution of the density profile over 
Q < t < IOO.OOOta, where the different colours represent dif- 
ferent times (as denoted by the colour bar). It is clear that the 
density profile is modified over time, and at f = 100, OOOt^ the 
density aty = is approximately 95.5% of its equihbrium value. 
Figure 5b shows the evolution of the temperature profile, and we 
see that there is a steady increase in temperature over the whole 
evolution, with maximum increase at y - 0. 
It is clear that the visco-resistive heating has substantially 
changed the background density profile, and that there is a sub- 
stantial flow of plasma in the positive )'-direction, causing a re- 
duction in the density since the boundary conditions do not allow 
for this plasma to be replaced by an in-flow. As can be seen from 
Figure 5a, it appears that the rate-of-change (i.e. the decrease) 
in the density is decreasing (this is confirmed below in §4.1). 
There are two possible explanations for this decrease in the rate- 
of-change: either the bulk flow (in the positive >'-direction) has 
moved so much density that an opposing pressure gradient (in 
the negative j'-direction) has been set up, or it could be that at 
later times there is less visco-resistive heating in our system. 
To examine the build-up of opposing pressure gradients, we can 
look at the evolution of dp/dy in our system (where p is total 
pressure), and this can be seen in Figure 5c. We can see that there 
are both positive and negative perturbations, related to the oscil- 
latory motions, but crucially we find no evidence for the build-up 
of an opposing pressure gradient in the negative y-direction. 
Let us now investigate the alternative hypothesis; at later times 
there is less visco-resistive heating in our system. Recall that 
equation (41) highlighted the role of kj in our system; kj is rel- 
evant both to the damping of the driven Alfven wave and to the 
enhanced bulk flow. Figure 5d shows the evolution of kj in our 
system, and we see that the magnitude of kj is decreasing over 
time. The behaviour of kj and the density profile are directly re- 
lated by equation (35). 

This provides the explanation for the decrease in the rate-of- 
change of the density profile: as the simulation proceeds, the 
visco-resistive damping heats the plasma, leading to a decrease 
in density. Thus (through equation 35), the magnitude of kj de- 
creases, thus there is less damping in the system, and thus the 
rate of visco-resistive heating decreases. This in turn leads to a 
decrease in the rate-of-change of density. There is a strong feed- 
back effect: a decrease in density leads to a decrease in \ki\ and 
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Fig. 5. Long-term evolution of (a) density, (b) temperature, (c) pressure gradient (dp/dy), and (d) k/, over < f < 100, OOOta, where the different 
colours represent different times. The colour bar denotes the different times. 
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Fig.6. (a) Transverse perturbations V; for 77 = v = 0.01, /3 = 0.1 plasma with po = 5, where red dashed lines represent the envelope V; = 
±Aexp{kiy). (h) Longitudinal perturbations Vy for same plasma, where black dashed line represents y = and blue horizontal line denotes 
\y = 0. In both subfigures, t = IOOOt^ and the vertical black dotted line represents y = v^f. 



thus a decrease in the nonideal heating, and thus a decrease in 
the bulk flow. 



4.1. Dependence of equilibrium density profile 

From equation (35) and the discussion in §4, it is clear that the 
larger the value of in our system, the more nonideal heating 



occurs at that time. Thus, there should be a strong dependence 
upon our choice of equilibrium density profile. 

Let us consider a ID system similar to that of §4 (i.e. /3 - 0.1, 
T] - V - 0.01, boundary conditions given by equation 31) where 
we set Po = 5. Figure 6a shows the behaviour of v~ in this system 
at f = lOOOr^i. As in Figure 2a, we see that the driven Alfven 
wave is damped. The magnitude of this resistive damping can 
be reproduced using dispersion relation (34), where we replace 
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Fig. 7. Long-term evolution of(a) density and (b) temperature, over < / < 100, OOOt^, where the dilferent colours represent different times. The 
colour bar denotes the different times. 
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Fig. 8. Evolution of dp/dt at^; = 1 for (a) po = 5 and (b) po = 1 plasmas. 
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Fig. 9. (a) Equilibrium density profile (black) and equilbrium Alfven speed profile (red), (b) Equilibrium temperature profile. 



77 by 7/ + V. Note that the wave is more rapidly damped than that 
in Figure 2a, since in the po = 5 system there is a larger value 
of \ki\ {\ki\ = 6.99 X 10"^^ as opposed to \ki\ = 6.25 x 10""* in the 
Po = 1 system). This can also be seen from equation (35). 

Figure 6b shows the behaviour of v,, at f = lOOOr^. As in Figure 
3a, we observe three components to the propagation: a pon- 
deromotive component, an acoustic (boundary-driven) compo- 



nent and a monotonically-increasing quadratic component, i.e. 
the bulk flow in the positive y-direction, as described previously. 
Note that the bulk flow is significantly stronger than that exam- 
ined in §4. This is because in the po = 5 system we have a larger 
value of \ki\ and thus more visco-resistive damping. Hence, there 
is more visco-resistive heating, which increases the local gas 
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pressure by a greater degree, and it is this thermal-pressure gra- 
dient that drives the bulk flow. 

Figures 7a and 7b show the long-term evolution of the density 
and temperature profile, respectively, in the po = 5 system. As 
in the po = 1 system, it is clear from Figure 7a that the density 
profile has decreased substantially during the simulation, and the 
density at y = is approximately 53.6% its initial value after 
t - 100, OOOt^. Figure 7a also clearly shows the bulk flow of 
density in the positive y-direction. Figure 7b shows that there is 
a monotonic increase in temperature over the whole simulation 
with maximum increase aty -0. 

Thus, it appears that a key variable in our system is the choice 
of equilibrium density profile: the greater the value of po, the 
greater the (initial) value of \kj\ and hence the greater the magni- 
tude of visco-resistive damping (which heats the plasma, which 
increases the thermal pressure, which drives the bulk flow due 
to pressure gradients). The bulk flow decreases the value of p, 
which decreases the value of \kj\, which reduces the amount 
of visco-resistive damping, and so on (strong feedback ef- 
fect). Thus, we can understand the nature of the bulk-flow phe- 
nomenon over a range of density values. 

Figures 8a and 8b show the rate-of-change of density (dp/dt) at 
y = 1 for the po = 5 and po = 1 systems, respectively. Let us first 
consider Figure 8a. it is clear that rate-of-change is decreasing as 
time evolves, i.e. \dp/dt\ is decreasing. Thus, as the simulation 
proceeds, density decreases at a slower rate, as expected from 
our explanation that as density decreases, there is less heating, 
leading to less thermal expansion. Figure 8b shows that, again, 
the rate-of-change of density is decreasing in the po = 1 system 
but at a much slower rate than that seen in Figure 8a. Again, this 
is in agreement with our explanation. 

5. Two-dimensional inhomogeneous plasma 

We now extend the model of §3 to include an inhomogeneous 
density profile. This will naturally introduce the phase mixing 
phenomenon into our system, in addition to the others akeady 
discussed. 

We consider the following equilibrium density profile; 



po(x) =1-1- 4sech" 



where, in this paper we set a - 35. This equilibrium density 
profile ( Po) can be seen in Figure 9a (black line); the red line 
represents the corresponding equilibrium Alfven-speed profile. 
The equilibrium temperature profile (Tq) is shown in Figure 9b, 
which has been chosen such that the equilibrium gas pressure is 
constant. 

As detailed in §2.3, we solve MHD equations (1) in a numer- 
ical domain of (x,y) e [-100, 100] x [0, 10"^] with a uniform / 
stretched grid in the x- / y-direction. We drive our system with 
linearly-polarised, sinusoidal Alfven waves at y = (boundary 
condition 31). All other quantities have zero gradient boundary 
conditions at the y-boundaries and we utilise periodic boundary 
conditions at the x-boundaries. 

The numerical results can be considered over two times scales; 
short-term evolution, over which the classical phase mixing so- 
lution dominates, and the long-term evolution over which non- 
linear efl'ects can become important. 

5. 1. Evolution ofw-: Alfven waves 

Figure 10 shows a contour plot of v- for < x < 100, < 
y < 500 at time f - IOOOta. This contour plot shows the prop- 
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Fig. 10. Contours of V; for < x < 100, < y < 500, at / = IOOOt^. 
White dashed line denotes x = 35. The temporal evolution is shown in 
the movie available in the on-line edition. 



agation of the Alfven waves in our system, where the light and 
dark bands denote the peaks and troughs, respectively, of the 
wave motions in the z-direction. Three different behaviours are 
apparent; the oscillations along x - correspond to the propa- 
gation of the Alfven wave in a po = 5 plasma and a cut along 
X = reveals a profile identical to Figure 6a. Similarly, the os- 
cillations along X - 100 correspond to the propagation of an 
Alfven wave in a po = 1 plasma and a cut along x = 100 is iden- 
tical to that in Figure 2a. As expected, the propagation along 
X = (where po[0,3'] = 5) is damped much more rapidly than 
that along x = 100 (where po[100,3'] = 1). Thus, the propaga- 
tion along X - and x - 100 and the (visco-resistive) damping 
mechanism can be fully understood from the corresponding one- 
dimensional results. 

Let us now consider a cut along x = 35 in Figure 10, i.e. 
Vj(35,)'); this can be seen in Figure 11. Here, we see that the 
wave is rapidly damped (more rapidly than standard visco- 
resistive damping). 

The classical phase mixing solution of Hey vaerts & Priest (1983) 
is given by; 



= +A exp (kjy) ■ exp 



4 









(43) 



where k\\ - Lo/vAix) and R = cij/{t] + v)^j^ (In^n)] . This solu- 
tion is overplotted in Figure 1 1 (red dashed line) and the agree- 
ment is excellent. Thus, at early times, the damping of V; is well 
understood as that corresponding to the classic phase mixing 
mechanism. 

A movie of the evolution of v, over the whole simulation (up 
to f = 100, OOOt/i) is associated with Figure 10. It is clear that 
the behaviour of v- changes very little over the whole simulation, 
and that Alfven waves are continuously driven into the numerical 
domain; there is no evidence of the system choking itself off. 
Note that the small changes in v. are dictated by the changes in 
the density profile (see §5.4 below). 
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Fig. 11. Plot of v.(35,;y). Heyvaerts & Priest solution overplotted 




Fig. 12. Shaded surface of Vy(x,y) at f = IOOOta. The temporal evolu- 
tion is shown in the movie available in the on-line edition. 



5.2. Evolution ofwyi Slow magnetoacoustic waves 

Figure 12 shows a shaded-surface of Vv(ji:,3') at f = IOOOta- As in 
the analysis of v , a cut along x - Q reveals behaviour identical 
to that of Figure 6b, and a cut along x = 100 reveals behaviour 
identical to that of Vy in a one-dimensional plasma with p - O.l, 
Tj = V - 0.01 at f = 1000t/i, i.e. similar to Figure 3a. Thus, as in 
our one-dimensional analyses, we see that our simulation con- 
tains three types of longitudinal motions: slow magnetoacoustic 
waves (previously called acoustic waves in §3), ponderomotive 
wave components (propagating at the Alfven speed) and a bulk 
flow in the positive y-direction. As expected, all amplitudes are 
of the order A^. 

A movie of the evolution of v,, over the whole simulation (up to 
t - 100, OOOta) is associated with Figure 12. It is clear that the 
growth of v^, saturates, and that the behaviour of v,, is dictated by 
that of V; and the heating profile. Thus, the complete evolution 
of Vy is of secondary importance in understanding the evolution 
of the density and temperature profiles. 

5.3. Evolution ofv,,: Fast magnetoacoustic waves 

In addition to Alfven waves and slow magnetoacoustic waves, 
there is now a third type of MHD wave present in our system: 
fast magnetoacoustic waves, which can be seen in Figure 13. 
These fast waves are nonlinearly generated by transverse gradi- 
ents in the Alfven-speed profile, in agreement with the analyt- 
ical work of Nakariakov et al. (1997). The fast waves are per- 
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Fig. 15. Density profiles of p(x, y) along the lines x = (blue), x = 35 
(black) and x = 100 (red), at t = 100, OOOr^. 



manently generated by phase mixing, and are refracted towards 
regions of lower Alfven speed. This can be clearly seen in Figure 
13a, which shows a contour plot of behaviour at f = IOOt^. 
These oscillations have a maximum around x - 35 and are re- 
fracted in towards x = (i.e. towards regions of lower Alfven 
speed, also see Figure 9a). These wave motions are generated 
with a frequency twice that of the driving frequency, which is in 
agreement with previous analytical predictions (Nakariakov et 
al. 1997; 1998). 

A shaded-surface of the behaviour of Vj^ at f = lOOOr^i can be 
seen in Figure 13b. Here, we can see that the generated fast 
waves have propagated across the magnetic fieldlines and have 
set-up an interference pattern. The interference pattern also re- 
sults from waves generated along x = 35 overlapping with those 
generated along x - -35, and there is also overlap due to the 
periodic boundary conditions. However, as first noted by Botha 
et al. (2000), these fast waves grow initially and then saturate. 
Botha et al. explain this saturation in terms of simple wave in- 
terference: physically, the saturation of fast waves is due to de- 
structive interference from incoherent sources. 
A movie of the evolution of v^ over the whole simulation (up to 
t - 100, OOOt/i) is associated with Figure 13. As with the evo- 
lution of Vy, it is clear that V y quickly saturates at amplitudes of 
order and, hence, it is only of secondary importance in under- 
standing the evolution of the density and temperature profiles. 

5.4. Evolution of density profile 

Let us now look at the evolution of the density profile in our sys- 
tem. Figure 14 compares shaded-surfaces of the initial density 
profile (Figure 14a) at f = and the final density profile (Figure 
14b) at f = 100, OQQta at the end of our simulation. Note that we 
have presented -100 < x < here, for clarity, and recall that the 
simulation is symmetric about x = 0. 

It is clear that the density profile has changed substantially dur- 
ing the simulation. Figure 15 shows cuts along x = 0, 35 and 
100 at the end of the simulation. The cuts along x = and 
X = 100 are in excellent agreement with those predicted by 
the one-dimensional model (i.e. compare to Figures 5a and 7a 
at f = 100, OOOt/i). Thus, the nature of the density change at 
these locations can be understood in terms of the bulk-flow phe- 
nomenon described in §4. 

The nature of the density change along x = 35 is however com- 
pletely different. Here, the bulk flow has been suppressed (very 





Fig. 14. (a) Shaded surface of p(x,y) at / = (i.e. equilibrium density profile), (b) Shaded surface of p(x,y) alt = 100, OOOta. 



little density change can be seen after y ^ 225) and instead the 
localised decrease in density is due to the strong, localised heat- 
ing from phase mixing (see below). 

5.5. Evolution of temperature profile 

Figure 16 shows contour plots of T - To at times (a) t = IQQQta, 
(b) t = 10,000ta and (c) t = 100,000ta. In Figure 16a, we 
see that there is strong, localised heating (maximum occurs at 
X - 38.6, y - 70.4) which results from the enhanced visco- 
resistive dissipation associated with phase mixing, producing a 
distinctive teardrop shape. Independently, the heating in the bot- 
tom left-corner is the visco-resistive heating associated Alfven- 
wave damping along pQ - 5. 

At a later time (Figure 16b, t - 10, OOOta) the overall tempera- 
ture profile appears qualitatively similar to that in Figure 16a, i.e. 
phase mixing is still the dominant heating mechanism. However, 
the maximum of T - To (located at x = 39.6, _v = 80.5) is 
now approximately ten times greater (note that the colour bar 
has changed). 



Figure 16c shows the temperature profile at the end of our simu- 
lation it - 100, OOOTyi). Again, qualitatively the contours appears 
similar to those in (a) and ib). However, the magnitude has again 
increased substantially (the colour bar has again changed) and 
the maximum of T - Tq is now located at x = 37.7, y - 108.8. 
Again, phase mixing is still the dominant heating mechanism. 
Figure 16d shows the location of the maximum of T - To for 
2000ta < t < 100,000ta (note the different length scales on 
the axes). It is clear that the location of the maximum of T - Tq 
changes during our simulation, and that this drifting is in the di- 
rection of decreasing Alfven-speed profile. This result addresses 
one of the fundamental question asked in §1, i.e. we find that 
drifting of the heating layer in nonideal, nonlinear simulations 
of phase mixing is possible. However, we have also shown that 
such drifting is very small and takes a very long time, at least for 
the parameters considered in this paper. 

Finally, note that before t - 2000t, the maximum of T - Tq 
sometimes occurs due to the bulk-flow heating, rather than due 
to phase mixing (phase mixing takes a finite amount of time to 
become the dominant heating mechanism) and thus, for clarity, 
we have not included these points in Figure 16d. However, for 
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comparison we have included the location of the maximum of 
T -Toatt = lOOOr (i.e. x = 38.6, y = 70.4) and this is denoted 
by a star. 

6. Conclusions 

We have investigated the nonlinear, nonideal behaviour of 
Alfven wave propagation and phase mixing within an inhomoge- 
neous environment, over long timescales. The governing MHD 
equations have been solved in ID and 2D environments using 
both analytical techniques and numerical simulations. 
In an ideal, one-dimensional study (no density inhomogeneity, 
d/dx = 0) we find that by driving a linear Alfven wave (in v,) 
into our numerical domain, we also generate two types of lon- 
gitudinal wave: boundary-driven acoustic waves (propagating at 
speed Cj) and a nonlinear perturbation (propagating at v^) which 
is driven by the ponderomotive force. Both these motions have 
a frequency twice that of the driven Alfven wave, and an ex- 
act mathematical solution for both wave types was derived. The 
acoustic wave is degenerate under the j8 = approximation (see 
Appendix A) but the ponderomotive wave is always present in a 
nonlinear system. 

We find that the addition of resistive and viscous effects nat- 
urally leads to visco-resistive damping (through the dispersion 
relation) and thus to visco-resistive heating, which in turn leads 
to the introduction of a new phenomenon: a bulk flow in the pos- 
itive >'-direction. Physically, this bulk flow is a direct response 
to the visco-resistive heating in the system: the heating increases 
the local temperature which increases the local thermal pressure. 
The resulting presure gradient drives the bulk flow. The bulk flow 
is purely in the positive y-direction because the boundary con- 
ditions are held fixed (v^ = 0). As a result of this out-flow, a 
significant reduction in density occurs since the boundary con- 
ditions do not allow for the plasma to be replaced by an in-flow. 
Note that the bulk flow is not due to the ponderomotive force, 
otherwise such a flow would be apparent in the ideal, nonlinear 
system. 

We find that the magnitude of the visco-resistive heating is 
strongly dependent upon kj and thus on our choice of po; the 

greater the value of po, the more visco-resistive heating that oc- 
curs. More visco-resistive heating leads to a stronger bulk flow 
in the longitudinal direction. 

We also investigated the effect of driving an Alfven wave in a 
nonideal plasma over very long timescales. We find that the equi- 
librium density profile changes substantially over the duration of 
our simulation (end of the simulation at / = 100, OOOta) specifi- 
caUy due to this bulk-flow phenomenon (which is itself directly 
due to the visco-resistive heating). For a po = 1 plasma, we find 
a decrease in density of about 4.6% at y = but for a po = 5 
plasma we find a decrease of about 46.4%. We also find that 
the rate-of-change of density at j = is decreasing as the sim- 
ulation proceeds. This is because as the density decreases, the 
value of \kj\ also decreases. This smaller value of \k;\ reduces 
the amount of visco-resistive damping, which slows the rate-of- 
change of density, and so on. Thus, there is a strong feedback 
effect in our system, and we can explain the nature of the bulk 
flow and change of density profile over a range of density val- 
ues. We conclude that the bulk-flow phenomenon is a natural 
consequence of driving an Alfven wave in a nonideal plasma. 
We then extended our analysis to include an inhomogeneous 
density profile in a two-dimensional system and, as before, we 
drive our system with sinusodial, linear Alfven waves. As ex- 
pected from Heyvaerts & Priest (1983), our simulation displays 
classic phase mixing. A cut showing Vj as a function of height 



along X = 35 provided an excellent match with the analytical 
predictions of Heyvaerts & Priest (1983). Cuts along x = Q 
( Po = 5) and x = 100 ( po = 1) were identical to those of 
the corresponding visco-resistive damped Alfven waves seen in 
our one-dimensional analysis. 

Our one-dimensional analysis also explained the longitudinal 
motions (Vy) in our system and we again observed the boundary- 
driven acoustic waves (interpreted as slow magnetoacoustic 
waves in 2D, with speed c,), our ponderomotive waves (prop- 
agating at va) and the bulk-flow phenomenon. 
In addition to Alfven waves and slow magnetoacoustic waves, 
our 2D system also contains fast magnetoacoustic waves (seen 
primarily in V;c). These fast waves are continuously generated by 
Alfven-wave phase mixing at a frequency twice that of the driven 
Alfven wave, and propagate across magnetic fieldlines and away 
from the phase mixing layer (as predicted by Nakariakov et al. 
1997; 1998). These waves are also refracted towards regions of 
lower Alfven speed. 

Since there is a permanent leakage of energy away from the 
phase mixing layer, it is possible that these fast waves can cause 
indirect heating of the plasma as they propagate away and dis- 
sipate far from the phase mixing layer itself, thus spreading the 
heating across the domain (Nakariakov et al. 1997). However, 
due to our choice of amplitude (A = 0.01) we find that only a 
smaU fraction of the Alfven-wave energy is converted into fast 
waves and, thus, only a small amount of indirect heating occurs. 
Over long timescales, we find that both the slow and fast waves 
saturate, and always remain at amphtudes of the order A^. 
Hence, these waves only play a secondary role in the heating 
of the plasma. 

We find that at the end of our simulation, the density profile has 
changed substantially. This is due to two effects: firstly, the bulk- 
flow phenomenon is naturally present in our system, with the 
density changes as predicted by our one-dimensional analysis. 
Secondly, there is a substantial decrease in density at the loca- 
tions of phase mixing heating. 

The change in the background density profile does not affect the 
propagation of the Affven wave (ampUtude of the order A) to a 
great degree, and the behaviour of Vx (fast waves) and Vy (slow 
waves) do not have a strong feedback effect on the system (both 
amplitudes of the order A^) and only play secondary roles in the 
system (i.e. A = 0.01, so feedback effect is only of the order 
of 1%). However, these conclusions may be different for Alfven 
waves driven with a much larger amplitudes. 
We have also investigated how the temperature profile changes 
during our simulation. We find a substantial increase in localised 
temperature, generated by the phase mixing mechanism, over the 
duration of our simulation, and find that the maximum tempera- 
ture increases by a factor of about a hundred during the simula- 
tion. We also find that the location of the maximum temperature 
drifts during our simulation, thus providing an answer to one of 
the key questions asked in §1, i.e. we find that drifting of the 
heating layer is possible in a nonlinear phase mixing scenario. 
However, we also find that this drifting is very small and occurs 
over a very long timescale (our simulation runs over 100, OOOta), 
at least for the parameters we have considered in this paper. For 
typical coronal parameters, ta = 60 seconds and so we are con- 
sidering timescales of approximately 69.4 days, which is too 
long to be physically important in the corona (i.e. other coronal 
processes act over much shorter timescales). However, it may 
be possible to increase the speed and magnitude of this drifting 
by significantly increasing the magnitude of heating in our sys- 
tem, i.e. attempt to increase the amount of enhanced dissipation 
from phase mixing. Thus, future studies could consider increas- 
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Fig. 16. Contour of T - To at (a) t = lOOOr^ ,(b)t= 10, OOOt^ and (c) f = 100, OOOta . Dashed white line denotes x = 35. (d) Plots (diamonds) the 
location of the maximum of T - To over 2000ta < t < 100, OOOta. Also plotted (single star) is the location at f = IOOOta. In subfigures (a) - (c), 
each colour bar represents magnitude, whereas in (d) the corresponding colour bar donotes time. 



ing the amplitude and frequency of the driven Alfven wave, or 
increasing the steepness of the gradient in our background den- 
sity inhomogeneity. 
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Appendix A: Alfven wave propagation in a 
one-dimensional, ideal, p = plasma 

Consider a cold (J3 - 0), ideal (77 = v = 0) plasma in our 
one-dimensional system (§3). Driving this system with bound- 
ary condition (31) generates an Alfven wave at y = 0; this can 
be seen in Figure A. la. The Alfven wave (v,) propagates in the 
direction of increasing/positive y with amplitude A - 0.01 and 
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there is no damping. The time of the snapshot can be read di- 
rectly from the y-axis using the relation t = y/vA- The ramp-up 
over the first four wavelengths is clearly visible. 
Figure A. lb shows the longitudinal motions (Vy) in the sys- 
tem. Here, v^, is driven by the nonlinear terms of equation (29) 
(~ B^dB^jdy) and again there is no damping. In this paper, we 
call this second-order nonlinear effect the ponderomotive ef- 
fect, and thus these longitudinal motions, which propagate at the 
Alfven speed, are driven by the ponderomotive force (Alfven 
wave-pressure gradients). At early times, these longitudinal mo- 
tions are governed by a simpUfied version of equation (29): 



1 d 
fipo dt 



(A.1) 



Assuming boundary conditions (31) in our ideal, cold plasma 
and that v^, v,. are initially zero, equation (A.l) has an analytical 
solution of the form: 



4v. 



{1 -cos(2w[?-);/va])} , 



(A.2) 



which is valid for < >■ < v^f. Note that these perturbations have 
an amplitude of 0{A^) and are always positive for an Alfven 
wave propagating in the positive y-direction (due to the con- 
stant of integration, and thus due to the boundary conditions). 
Finally, note that the frequency of these longitudinal motions is 
twice that of the driving frequency, and that the motions do not 
grow with time. Interestingly, the ramp-up to maximum ampli- 
tude now occurs over eight wavelengths, i.e. twice that of the 
driven wave. 



and the second has the form: 

V 2 , ,2 d2 4 2 



[dyj - 4vl J 
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(B.4) 



The form of equations (B.2) - (B.4) lead to equations (38) - (41). 
Equations (B.2) - (B.4) can also be derived explicitly from our 
choice of Vj-, i.e.: 

v^ = A sin {cot - kRy)e'''^ 

where as before k = kR + ikj. Thus (from equation 7) fo^ has the 
form: 



loBqA 



[ku sin (a)t - kny) 



yl{kl + k^) 

-k, cos (cot - kR}')] e*'^ . (B.5) 
Now we can explicitly determine the two contributions: 



dy 
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"A 
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y\{kl + k]) 
+kis\nl{ojt-kRy)\ e^'^ 
and the second has the form: 



[kR cos 2 (cjt - kRy) 



(B.6) 



Appendix B: Derivation of equations 38 - 41 

Equation (29) governs the driven, longitudinal motions in the 
system, and in ID {d/dx = 0) has the form: 



1 d 



B. 



dB, 



(r-i)77 d 



dB, 

dy 



(B.l) 



jupo dt\ dy j hpq dy 
Let us assume that our linear Alfven wave can be represented as: 
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(B.7) 



Equations (B.5) - (B.7) give equivalent solutions to equations 
(B.2 - B.4), and thus also lead to equations (38) - (41). 
Note that in both derivations, we have assumed an infinite har- 
monic solution for v^, whereas our numerical simulations are for 
driven harmonic wavetrains. Thus, our solutions are only vaUd 
for Cst <y < VAt. 



where A is the amplitude of our wave, co is the driving frequency, 
k - kR + ikj is our wavenumber, with complex conjugate k* = 
kR - ikj. Thus, from equation (7) (with A^3 = 0) we have: 



cdBqA 



]_^i(ul-ky) |_g-(M-ry) 

ik ik* 



(B.2) 



The right-hand-side of equation (B.l) has two contributions, 
which can be calculated using our expression for b,. The first 
term has the form: 
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Fig. A.l. (a) Transverse perturbations Vj for ideal = v = 0) j8 = plasma, {b) Longitudinal perturbations V3, for ideal /3 = plasma. In both 
subfigures, t = IOOOta and dotted line represents y = vaI. 



